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Abstract An expression for the stress tensor near an exter- 
nal boundary of a discrete mechanical system is derived ex- 
plicitly in terms of the constituents' degrees of freedom and 
interaction forces. Starting point is the exact and general 
coarse graining formulation presented by Goldhirsch in [I. 
Goldhirsch, Gran. Mat., 12(3):239-252, 2010], which is con- 
sistent with the continuum equations everywhere but does 
not account for boundaries. Our extension accounts for the 
boundary interaction forces in a self-consistent way and thus 
allows the construction of continuous stress fields that obey 
the macroscopic conservation laws even within one coarse- 
graining width of the boundary. 

The resolution and shape of the coarse-graining function 
used in the formulation can be chosen freely, such that both 
microscopic and macroscopic effects can be studied. The 
method does not require temporal averaging and thus can 
be used to investigate time-dependent flows as well as static 
and steady situations. Finally, the fore-mentioned continu- 
ous field can be used to define 'fuzzy' (highly rough) bound- 
aries. Two discrete particle method (DPM) simulations are 
presented in which the novel boundary treatment is exem- 
plified, including a chute flow over a base with roughness 
greater than a particle diameter. 

Keywords Coarse graining • Averaging ■ Boundary 
treatment ■ DPM (DEM) ■ Discrete mechanical systems ■ 
Homogenisation ■ Stress • Continuum mechanics • Granular 
systems 
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1 Introduction 

The main topic of this paper is the issue of coarse-graining, 
near a boundary. We consider the bulk method described by 
Isaac Goldhirsch 1 1 1, and extend it to account for boundary 
forces due to the presence of a wall or base. The Goldhirsch 
special edition of Granular Matter is an appropriate place 
to present some of the ideas that we have developed in this 
area. 

Continuum fields often need to be constructed from from 
discrete particle data. In molecular dynamics (T] and gran- 
ular systems [3^, these discrete data are the positions, ve- 
locities and forces of each atom or particle. In contrast, in 
the case of smooth particle hydrodynamics 15], the contin- 
uum system itself is approximated by a discrete set of fluid 
parcels. In all these methods, a crucially important issue is 
how to compute the continuum fields in the most appropriate 
way. Several techniques have been developed to calculate 
the continuum fields, see IfTOll and references therein. Partic- 
ularly the stress tensor is of interest: the techniques include 
the Irvin-Kirkwood's approach f6l or the method of planes 
|7j. Here, we use the coarse-graining approach (CG) as first 
described in JSl. 

The CG method |[8][I] has several advantages over other 
methods, including: i) the fields automatically satisfy the 
conservation equations of continuum mechanics; ii) it is not 
assumed that the particles are rigid or spherical, and Hi) the 
results are valid for single particles (no averaging over en- 
sembles of particles is required). The only assumptions are: 
each particle pair has a single point of contact, the contact 
area can be replaced by a contact point, and collisions are 
not instantaneous. 

In Sect. 121 we use the derivation of |[T1 to extend the CG 
method to account for the presence of a boundary. Explicit 
expressions for the resulting continuum fields are derived. In 
Sect. 12.51 an alternative stress definition is proposed extend- 
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ing the stress field into the boundary region. In Sect. |3] the 
approach is tested with two DPM simulations, and in Sect. 
|4]we draw conclusions. 



2 Theory 

2.1 Assumptions and notation 

We are interested in deriving macroscopic fields, such as 
density, velocity and the stress tensor from averages of mi- 
croscopic variables such as the positions, velocities and for- 
ces of the constituents. Averaging will be done such that the 
continuum fields, by construction, satisfy conservation laws. 
Vectorial and tensorial components are denoted by Greek 
letters in order to distinguish them from the Latin particle- 
indices Bold vector notation will be used when appro- 
priate. We will follow the derivation of fT), but extend it 
by introducing two types of particle: flowing particles 
{1,2,...,A^} and boundary particles {N + I, . . . ,N + K}. 

Each particle / has mass m,-, center of mass position ria 
and velocity v,a . The force /,« acting on particle ; is a com- 
bination of the sum of the interaction force fija with another 
particle j, the interaction force with a boundary particle 
k, and a body force bja (e.g., gravity). 



N 

fia = 52 -^V" 



N+K 

V, fika 

k=N+\ 



-b,a, i<N. 



(1) 



The interaction forces are binary and anti-symmetric such 
that action equals reaction, fija — —fjta, i-J < N. We as- 
sume that each particle pair i < N, j < N + K has, at 
most, a single contact point, Cija, at which the contact forces 
act. The positions of the boundary particles are fixed, as if 
they had infinite mass. The trajectories of the flowing parti- 
cles are governed by Newton's second law and if tangential 
forces and torques are present, rotations follow from the an- 
gular form of Newton's law. 

In the following sections, we commence from Ref. 1 1 1 
to derive definitions of the continuum fields. To be precise, 
a body force density is introduced to account for body forces, 
and to incorporate boundary effects an interactionforce den- 
sity (IFD) is introduced. While the idea of an IFD is more 
generally applicable (e.g., for mixtures), it is employed here 
to account for the presence of a boundary. 



2.2 Coarse graining 

From statistical mechanics, the microscopic mass density of 
the flow at a point at time t is defined by 



p'^'%r,t) = Y^m,8{r-r,{t)), 

1=1 



(2) 



where 5{r) is the Dirac delta function. We use the following 
definition of the macroscopic density. 



p{r,t) = Y.^iW{r-ri{t)), 

i=l 



(3) 



i.e., we have replaced the Dirac delta function by an inte- 
grable 'coarse-graining' function W whose integral over the 
domain is unity. 

2.3 Mass balance 

The coarse-grained momentum density is defined by 



Pa{r,t) = Y^miViaW{r-ri). 

1=1 



(4) 



Hence, the macroscopic velocity field Va{r,t) is defined as 
the ratio of momentum and density fields, Va{r,t) = 
Pa{r,t)/ p{r,t). It is straightforward to confirm that pa and 
Pa satisfy the continuity equation (c.f. |[l]|8]), 

dp , dpa 



dt dra 



= 0. 



(5) 



2.4 Momentum balance 

Subsequently, we will consider the momentum conservation 
equation with the aim of establishing the macroscopic stress 
field. Gap ■ As we want to describe boundary stresses as well 
as internal stresses, the boundary interaction force density 
(IFD), ta, has been included, as well as the body force den- 
sity, ba, which are not present in the original derivation, IT]. 
The desired momentum balance equations take the form. 



dpa 



da, 



a/3 



dt drp drp 



-ta + bo 



(6) 



To determine the stress it is required to compute the tempo- 
ral derivative of (HJl, 

-Y.fiaW{r-r,) + Y^miVia^W{r~r,), (7) 



dt 



i=\ 1=1 

where /,« — m,dv,Q;/df is the total force on particle /. Using 
([T]), the first term in ^ can be expanded as 

N N N N+K N 

A„ = ^ E /'-/■«^ + L L + (8) 

i=\j=\J^i (=1 k=N+ 1 (=1 

with the abbreviation Wi = W{r — r,). The first term, which 
represents the bulk particle interactions, satisfies 

N N N N 

L L fija^ = I L fj,a^j 

N N 

= L L f'j-^j^ (9) 
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since fija — ^f jia and because the dummy summation in- 
dices can be interchanged. It follows from (|9]) that 

N N ^ N N 

i=lj=i+l 

Substituting ([TOll into dSJ yields 



(10) 



Expression (fTST i can be satisfied by defining the last term 
on the right hand side as the IFD. This however has the dis- 
advantage that the boundary IFD is located around the center 
of mass of the flowing particles. The more natural physical 
location of the boundary IFD would be at the interface be- 
tween the flowing and boundary particles. 

Therefore, we move the IFD to the contact points, c,<.a, 
between flowing and boundary particles: similar to ( fT2b . 



N N N N+K 

Aa = Y. L /'./«(^-^7) + L L fika^i + ba, dD 

/=l;=i+l i=U-=A'+l 

where ba = Jlibia^i is the body force density. 

Next, Aa is rewritten using Leibnitz's rule to obtain a 
formula for the stress tensor. The following identity holds 
for any continuously differentiable coarse-graining function W 

1 d 



(12) 



Jo OS 



{r-ri + srij)ds, 



where rija = ria — rja is the vector from rja to r,a. Substi- 
tuting identities ( fT2] i into (fTTT i yields 

d N N „i 



i=\j=i+\ 

N N+K 
i=lk=N+l 







(13) 



In Ref. m, it is shown that the second term in Q can be 
expressed as 



^ d 

E "^i^'ia 



dt 



_d_ 

dra 



(14) 



where v^^ is the fluctuation velocity of particle /, given by 
v\„ir,t)=Viait)-Vair,t). (15) 
Substituting ( fTST l and ( fT4b into momentum balance (|6]l yields 



drp 



drp 



''P i=lk=N+l 



where the kinetic and bulk contact contributions to the stress 
tensor are defined as 

N 

i=l 

<p = -tt fuanjp t y^{r - ri + srtj) ds. (17b) 

;= 0=1+1 ■'0 

Here, the underlined terms in (fT6l l are not in the original 
derivation presented in Ref. 1 1] and account for the presence 
of the boundary. 



ara Jo 



{r-ri + saik)ds, 



(18) 



where Wik = W{r - c,i) and = r,a - c,ia. Substituting 
( fTSl l into the last term in (fTSI l we obtain 

N N+K N N+K 

E E /;ta^ = E E f'ka^^'k 

i=\k=N+\ i=\k=N+\ 



_d_ 

dra 



N N+K 



i=ik=N+l 



E E fika^ikp / ^ {r - Ti + SOik) ds 



(19) 



Thus, substituting (fT9] l into ( fTSl l, we define the stress by 

(20a) 



where the contribution to the stress from the contacts be- 
tween flow and boundary particles is 

N N+K 



^ap = -. 

i=lk=N+l 

and the IFD is 

N N+K 

fa = E E fika'^ik- 

i=lk=N+l 



E E f'kaaikp / W{r-r,+sa,k)ds, (20b) 



(21) 



Equations (l20l l differ from the standard result of fT] by 
an extra term, a^'^, that accounts for the additional stress 
created by the interaction of the boundary with the flow. 
The definition ( 1211 1 gives the boundary IFD applied by the 
flowing particles; i.e., it has been constructed such that in 
the limit w — > 0, the IFD acts at the contact points between 
boundary and flow. 

Note, that this framework is general and can be used to 
compute more than boundary IFDs. For example, one can 
obtain the drag between two different species of interact- 
ing particles by replacing the flowing and boundary particles 
with the particles of the two species in the definition of the 
continuum fields. By placing an IFD at the contact points, 
the IFDs of both species are exactly antisymmetric and thus 
disappear in the momentum continuity equation of the com- 
bined system. In mixture theory, e.g. ifTTI . such interaction 
terms appear in the governing equations for the individual 
constituents and are called interaction body forces. These in- 
teraction body forces are an exact analog to the IFDs. There- 
fore, our approach can interpreted as treating the system as 
a mixture of boundary and flow particles and the IFD is the 
interaction body force between different species of particle. 
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Further, we note that the integral of the stress in ( [TtT i 
and ( |20| ) over the domain Q satisfies the virial definition of 
mechanical stress, 

Jn 

N N N N+K 

"1111 f'Janjp-Y, H fikaOikp- (22) 

i=l;=!+l 1= 1^=^+1 



2.5 Extending the stress profile into a base or wall 

In contrast to the previous subsection, an alternative stress 
definition is presented here, where the IFD and the stress are 
combined into a single tensor. Similar to ( fTSl i and ( fTSb . the 
following identity holds, 

Jq OS 

^nkp^ ^{r-n + sni,)ds, (23) 

orp Jo 

since the coarse-graining function W satisfies W {\r \ ^ °o) = 
0. Substituting (l23l l into ( fT6b we can obtain an alternative 
solution with zero IFD, = 0, where the stress is given by 

<p' = - L iT f.kankli r^ir- n + sr,k) ds. (24) 

This stress definition is not identical to the one in ( |20] | and 
( flTl l. It eliminates the IFD term entirely and provides a nat- 
ural extension of the stress into the boundary. However, the 
extended stress does contain contributions from both inter- 
nal and external forces, and the spatial integral of the stress 
components has to be extended to infinity. In singular spe- 
cial cases this can lead to artificial results. Another disad- 
vantage of Eq. ( l24b is its difficult interpretation due to the 
long-ranging integral. One could see it as the stress inside a 
'virtual/fake' wall-material on which the body-force is not 
acting (equivalent to foam with zero mass-density). How- 
ever, this is far fetched and not realistic, so that we rather 
stick to the formulation in Sect. 12.41 

It is also possible to extend the stress tensor to the bound- 
ary by other means, such as mirroring the stress at the bound- 
ary, or using a one-sided coarse-graining function. This is 
not discussed further since the first method requires a defi- 
nition of the exact location of the boundary, while the second 
method can introduce a spatial shift in the stress field due to 
spatial inhomogeneities. 



3 Results 

3.1 Contact model 

For illustrational purposes, we simulate a granular system 
with contact interaction forces. The statistical method, how- 
ever, is based only on the assumptions in fj2] and therefore 
can be applied more generally. We use a viscoelastic force 
model with sliding friction as described in detail in I^H. 
The parameters of the system are nondimensionalised such 
that the flow particle diameters are = 1, their mass m = 1, 
and the magnitude of gravity g—l - The normal spring and 
damping constants are A:" = 2 • 10^ and y" = 50, respectively; 
thus, the collision time is tc — 0.005 yjdjg and the coeffi- 
cient of restitution is e = 0.88. The tangential spring and 
damping constants are W = 2 /Ik" and y' = y", such that 
the frequency of normal and tangential contact oscillation, 
and the normal and tangential dissipation are equal. The mi- 
croscopic friction coefficient is set to ji^ ~ 0.5. We inte- 
grate the resulting force and torque relations in time using 
the Velocity- Verlet algorithm with a time step At — ft/50. 

We take the coarse-graining function to be a Gaussian of 
width, or variance, w. Other coarse-graining functions are 
allowed, but the Gaussian has the advantage that it produces 
smooth fields and the required integrals can be performed 
exactly. 

3.2 Quasi-static example in two dimensions 

In order to visualise definitions (l20l l and (1211 1. we firstly 
consider a two-dimensional configuration consisting of five 
fixed boundary particles and five flowing bulk particles, with 
gravity in the z-direction, see Fig.[T] The flow is relaxed un- 
til the flowing particles are static; hence, the only contribu- 
tion to the stress is due to the enduring contacts. To visu- 
alise the spatial distribution of the IFD and stress, the norms 
\t\ = and \(j\ — max|;[.|=i \cx\ (the maximum absolute 
eigenvalue), are displayed in Fig. [1] A very small coarse 
graining width, w = ii/8, is chosen to make the spatial aver- 
aging visible: the IFD, Eq. i2T[ . centers around the contact 
points between flowing and static particles, r,,ta, while the 
stress, Eqs. (|20| ) and (fTTI l. is distributed along the contact 
lines, narja and naCita- 

3.3 Three-dimensional steady chute flow 

Secondly, we consider a three-dimensional simulation of a 
steady uniform granular chute flow, see Fig.|2]and Ref. lfT2ll . 
The chute is periodic in the x- and y-directions and has di- 
mensions G [0,20] X [0,10]. The chute is inclined at 
Q ~ 26° and the bed consists of a disordered, irregular bound- 
ary created from fixed particles with size dbase = 2. The 
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Fig. 1 Grey circles denote a two-dimensional configuration of free and fixed boundary particles, with gravity in z-direction. Fixed particles are 
marked with a cross in the center. Contour plots show the spatial distribution of the norms of the boundary IFD and the contact stress (magnitude 
of largest eigenvalue). A very small coarse graining width, w = is chosen to make the spatial averaging visible: the IFD centers around the 
contact point, while the stress is distributed along the contact lines. 




Fig. 2 Steady chute flow over a very rough frictional surface of incli- 
nation 6 = 26° for A' = 1000 flowing particles. Gravity direction g and 
coordinates {x,y,z) as indicated. The domain is periodic in the x- and 
y-directions. Shade indicates speed; dark is slow and bright is fast. 



Fig. 3 Downward nornial stress ov, without (dashed) and with (solid) 
correction by the boundary IFD for w = d/4. The stress and IFD ex- 
actly match the weight of the flow above height z (red dotted), as ex- 
pected for steady flows. Grey lines indicate bed and surface location. 



chute contains 1000 flowing particles, which are initially 
randomly distributed. The simulation is computed until a 
steady state is reached. A screen shot of the steady-state sys- 
tem is given in Fig. |2] 

Depth profiles for steady uniform flow are obtained by 
averaging with a coarse-graining width w = d/4 over x E 
[0,20],ye [0,10] andfe [2000,2100]. The spatial averaging 
is done analytically, while we average in time with snapshots 
taken every tc/2. 

Note that the stress definitions (|20|) and (|24| satisfy 

^ = ^+.„. (25) 
After averaging in x, y and / directions, this yields 

~5 — = ^ 1"'"- *-26) 

ar^ or. 

Integrating over {z,°°), we obtain 

(^'az = <^az- / tadr^. (27) 



Thus, setting a = z in (l27b . the extended stress component 
a'., can be obtained without computing the semi-infinite line 
integral. 

The depth profile for the downward normal stress a^^ is 
shown in Fig. [3] Since the rough boundary is not at a fixed 
height, the stress gradually decreases at the bottom due to 
the decreasing number of bulk particles near the base. Due 
to the coarse graining, the stress tensor has a gradient even in 
the case of a flat wall, but the gradient disappears as w — > 0. 
Using the extended stress definition, the bed and surface lo- 
cations can be defined as the line where the downwards nor- 
mal stress a'-- vanishes and where it reaches its maximum 
value (to within 2%), see Fig. [3] Additionally, since the flow 
is steady and uniform, (|6]l yields the so called lithostatic bal- 
ance, a', — — /~ gzP dr^, which is satisfied with good accu- 
racy. 

4 Conclusions 

We have derived explicit expressions for the stress tensor 
and the interaction force density (IFD) near an external boun- 
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dary of a discrete mechanical system. These expressions were 
obtained by coarse-graining the microscopic equations and 
therefore exactly satisfy the governing balance laws of mass 
(O and momentum A boundary IFD was computed us- 
ing the contact points between the flow and the basal par- 
ticles. Our results can be extended to other IFDs, for ex- 
ample, the drag between two different species of particles. 
The power of our extension to Goldhirsch's method has been 
demonstrated by computing stress profiles for a chute flow 
over a fuzzy boundary. It avoids the problems inherent in 
other methods and gives the expected linear Uthostatic pro- 
file all the way to the base. 

The present formulation for boundary interaction forces 
allows us to draw the analogy to electrostatics, where the 
divergence of the electric field (analogous to the divergence 
of stress) is compensated by a charge-density source like our 
interaction force density ( |2T] ). The analogy can also be made 
to mixture theory where, by treating the system as a mixture 
of boundary and flow particles, the IFD is then interpreted 
as the interaction body force between the two species. 
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